###############
### Setup R ###
###############

### Clear terminal
cat("\014")

### Clear space
rm(list = ls())

## This is a slightly modified version of Michal Kolesar's BM_StandardErrors.R
## The original file was downloaded from:
##   https://github.com/kolesarm/Robust-Small-Sample-Standard-Errors
##   (Latest commit 2a80873 on Aug 26, 2015)
## We made minor changes for efficiency, as well as two changes that affect the
## code's functionality:
## 1) MatSqrtInverse() now stops with an error message if its argument is not of
##    full rank:
##    "Bell-McCaffrey SE undefined. This happens, e.g., when a dummy regressor is 1
##     for one cluster and 0 otherwise."
## 2) The list returned by BMlmSE() now includes "se".

# The MIT License (MIT)
#
# Copyright (c) 2015 Michal Kolesar
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

## Compute Bell-McCaffrey Standard Errors
library(sandwich)
library(Matrix)

message1 <- paste0(
  'Bell-McCaffrey SE undefined. This happens, e.g., when a dummy regressor is 1 ',
  'for one cluster and 0 otherwise.'
)

MatSqrtInverse <- function(A) {
  ##  Compute the inverse square root of a matrix
  if (rankMatrix(A) < NROW(A)) stop(message1)
  ei <- eigen(A, symmetric = TRUE)
  d2 <- 1/sqrt(ei$values)
  ## diag(d2) is d2 x d2 identity if d2 is scalar, instead we want 1x1 matrix
  ei$vectors %*% (if (length(d2)==1) d2 else diag(d2)) %*% t(ei$vectors)
}

BMlmSE <- function(model, clustervar=NULL, ell=NULL, IK=TRUE) {
  X <- model.matrix(model)
  sum.model <- summary.lm(model)
  n <- sum(sum.model$df[1:2])
  K <- model$rank
  XXinv <- sum.model$cov.unscaled # XX^{-1}
  u <- residuals(model)

  df <- function(GG) {                # Compute DoF given G'*Omega*G
    sum(diag(GG))^2 / sum(GG * GG)
  }

  if(is.null(clustervar)) {           # no clustering
    Vhat <- vcovHC(model, type="HC2")
    Vhat.Stata <- Vhat*NA

    M <- diag(n)-X %*% XXinv %*% t(X)       # annihilator matrix
    GOG <- function(ell) {           # G'*Omega*G
      Xtilde <- drop(X %*% XXinv %*% ell / sqrt(diag(M)))
      crossprod(M * Xtilde)
    }
  } else {
    if(!is.factor(clustervar)) stop("'clustervar' must be a factor")

    ## Stata
    S <- length(levels(clustervar)) # number clusters
    uj <- apply(u*X, 2, function(x) tapply(x, clustervar, sum))
    Vhat.Stata <- S/(S-1) * (n-1)/(n-K) * sandwich(model, meat = crossprod(uj)/n)

    ## LZ2
    tXs <- function(s) {
      Xs <- X[clustervar==s, , drop=FALSE]
      MatSqrtInverse(diag(NROW(Xs))-Xs%*% XXinv %*% t(Xs)) %*% Xs
    }
    tX <- lapply(levels(clustervar), tXs) # list of matrices

    tu <- split(u, clustervar)
    tutX <- sapply(seq_along(tu),function(i) crossprod(tu[[i]],tX[[i]]))
    Vhat <- sandwich(model, meat = tcrossprod(tutX)/n)

    ## DOF adjustment
    tHs <- function(s) {
      Xs <- X[clustervar==s, , drop=FALSE]
      index <- which(clustervar==s)
      ss <- outer(rep(0,n),index)     # n x ns matrix of 0
      ss[cbind(index,1:length(index))] <- 1
      ss-X %*% XXinv %*% t(Xs)
    }
    tH <- lapply(levels(clustervar), tHs) # list of matrices

    Moulton <- function() {
      ## Moulton estimates
      ns <- tapply(u, clustervar,length)
      ssr <- sum(u^2)
      rho <- max((sum(sapply(seq_along(tu), function(i)
        sum(tu[[i]] %o% tu[[i]])))-ssr) / (sum(ns^2)-n), 0)
      c(sig.eps=max(ssr/n - rho, 0), rho=rho)
    }

    GOG <- function(ell) {
      G <- sapply(seq_along(tX),
                  function(i)  tH[[i]] %*% tX[[i]] %*% XXinv %*% ell)
      GG <- crossprod(G)

      if (IK==TRUE) {            # IK method
        Gsums <- apply(G, 2, function(x) tapply(x, clustervar, sum)) # Z'*G
        GG <- Moulton()[1]*GG + Moulton()[2]*crossprod(Gsums)
      }
      GG
    }
  }

  if (!is.null(ell)) {
    se <- drop(sqrt(crossprod(ell,Vhat) %*% ell))
    dof <- df(GOG(ell))
    se.Stata <- drop(sqrt(crossprod(ell,Vhat.Stata) %*% ell))
  } else {
    se <- sqrt(diag(Vhat))
    dof <- sapply(seq(K), function(k) df(GOG(diag(K)[,k])))
    se.Stata <- sqrt(diag(Vhat.Stata))
  }
  names(dof) <- names(se)

  return(list(vcov=Vhat, dof=dof, adj.se=se*qt(0.975,df=dof)/qnorm(0.975),
              se=se,
              se.Stata=se.Stata))
}


### Load packages
library(rio)
library(nnet)
library(lmPerm)

### Print system info
sessionInfo()

###########################
### Setup Missions Data ###
###########################

### Import data
dat.sub <- rio::import("tcef.csv")
colnames(dat.sub) # Check column names

### Check for balance using multinomial logit, as reported on page 13
#### Estimate null model
multi.out.0 <- multinom(tr ~ 1, data = dat.sub)
#### Save LL
LL.r <- logLik(multi.out.0)[1]
#### Estimate full model
multi.out.1 <- multinom(tr ~ as.factor(ccode_host), data = dat.sub)
summary(multi.out.1)
#### Save LL
LL.u <- logLik(multi.out.1)[1]
#### Calculate DF
degfree <-
  length(summary(multi.out.1)$standard.errors) -
  length(summary(multi.out.0)$standard.errors)
##### LR test statistic and p-value
2*(LL.u - LL.r)
1 - pchisq((2*(LL.u - LL.r)), df = degfree) # p

### Calculate reply rates, as reported on page 14
round(mean(dat.sub$reply_prisoner)*100, 2) # 29.30%
round(mean(dat.sub$reply_baseline)*100, 2) # 24.36%

### Calculate reply rates conditional on treatment status, as reported in Table 2 on page 14
round(mean(dat.sub$reply_prisoner[dat.sub$tr==0]), 3) # Sanctioning = 0, Norms = 0; 0.247
round(mean(dat.sub$reply_prisoner[dat.sub$tr==1]), 3) # Sanctioning = 1, Norms = 0; 0.256
round(mean(dat.sub$reply_prisoner[dat.sub$tr==2]), 3) # Sanctioning = 0, Norms = 1; 0.46
round(mean(dat.sub$reply_prisoner[dat.sub$tr==3]), 3) # Sanctioning = 1, Norms = 1; 0.216
#### Compare control response rate to norm response rate, as reported on page 14
(round(mean(dat.sub$reply_prisoner[dat.sub$tr==2]), 3)/round(mean(dat.sub$reply_prisoner[dat.sub$tr==0]), 3))*100

### Conduct omnibus chi-sq test, as reported on page 14
chisq.test(xtabs(~dat.sub$reply_prisoner+dat.sub$tr))

### Estimate LPM for results reported in Figure 4 and in Appendix A in Table 3
#### Model includes treatment indicators, the blocking covariate, fixed effects for other randomized aspects of our email texts, and a count of mission types
##### Estimate model
ols.fit <-  lm(dat.sub$reply_prisoner ~ as.factor(dat.sub$sanctioning_no_norms) +
                 as.factor(dat.sub$norms_no_sanctioning) + as.factor(dat.sub$norms_sanctioning) +
                 as.factor(dat.sub$sender_two) + as.factor(dat.sub$wave) +
                 dat.sub$type + as.factor(dat.sub$ccode_host) + as.factor(dat.sub$separator) +
                 as.factor(dat.sub$separator_two) + as.factor(dat.sub$valediction), singular.ok = TRUE)
##### Create matrix to save estimates
estimates <- matrix(nrow = length(ols.fit$coefficients)-1, ncol = 6)
##### Create column with variable names
estimates[,1] <- as.vector(names(ols.fit$coefficients)[-1])
##### Compute BM SE and save them
bm <- BMlmSE(model = ols.fit)
##### Save LPM point estimates
point.estimate <- coef(ols.fit)
##### Set critical value
critical.value <- qt(0.975, df = bm$dof)
##### Set margin of error
margin.of.error <- critical.value * bm$se
##### Save CI
ci = c(point.estimate - margin.of.error, point.estimate + margin.of.error)
point.estimate  # Estimated average treatment effect
bm$se           # HC2 robust SE
bm$dof          # Bell-McCaffrey degrees of freedom
ci              # Bell-McCaffrey confidence interval
##### Save point estimates, SE, and t-stat
estimates <- cbind(point.estimate,bm$se,point.estimate/bm$se)[2:4,]
estimates # These results are reported in Appendix A in Table 3

### Create Figure 4
par(mfrow=c(1,1))
var.names <- c("Shaming", "Norms", "Norms and Shaming")
par(
  oma = c(0,0,0,0), # Since it is a single plot, I set the outer margins to zero.
  mar = c(5,4,5,4) # margins adjusted to reflect changing the locations of the labels
)
# create an empty plot for total customization
plot(NULL, # create empty plot
     ylim = c(-.3, .3), # set xlim by guessing
     xlim = c(1, .4+length(var.names)), # set ylim by the number of variables
     axes = F, xlab = NA, ylab = NA)
abline(h=0, lty=3)
est <- as.numeric(estimates[,1]) # conveniently store the estimates (minus the constant)
se <- as.numeric(estimates[,2]) # conveniently store the std. errors (minus the constant)
for (i in 1:length(est)) { # loop over a counter the length of the estimate vector
  lines(c(.2+i, .2+i), c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), lwd = 2, col="#CC6666")
  lines(c(.2+i, .2+i), c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), lwd = 4, col="#CC6666")
  points(.2+i, est[i], pch = 19, cex = 2, col="white")
  points(.2+i, est[i], pch = 19, cex = 1, col="black")
  text(.2+i, est[i], var.names[i], xpd = T, cex = .8, pos = 3, offset=5) # add variable labels above the points
}
axis(side = 2) # add bottom axis
mtext(side = 2, "Estimated Treatment Effects", line = 2.5)

### Estimate LPM for results reported in Figure 5 and in Appendix B in Table 4
#### Model includes treatment indicators, the blocking covariate, fixed effects for other randomized aspects of our email texts, and a count of mission types
##### Estimate model
ols.fit.helpful <-  lm(dat.sub$helpful_high ~ as.factor(dat.sub$sanctioning_no_norms) +
                 as.factor(dat.sub$norms_no_sanctioning) + as.factor(dat.sub$norms_sanctioning) +
                 as.factor(dat.sub$sender_two) + as.factor(dat.sub$wave) +
                 dat.sub$type + as.factor(dat.sub$ccode_host) + as.factor(dat.sub$separator) +
                 as.factor(dat.sub$separator_two) + as.factor(dat.sub$valediction), singular.ok = TRUE)
##### Create matrix to save estimates
estimates.helpful <- matrix(nrow = length(ols.fit.helpful$coefficients)-1, ncol = 6)
##### Create column with variable names
estimates.helpful[,1] <- as.vector(names(ols.fit.helpful$coefficients)[-1])
##### Compute BM SE and save them
bm <- BMlmSE(model = ols.fit.helpful)
##### Save LPM point estimates
point.estimate <- coef(ols.fit.helpful)
##### Set critical value
critical.value <- qt(0.975, df = bm$dof)
##### Set margin of error
margin.of.error <- critical.value * bm$se
##### Save CI
ci = c(point.estimate - margin.of.error, point.estimate + margin.of.error)
point.estimate  # Estimated average treatment effect
bm$se           # HC2 robust SE
bm$dof          # Bell-McCaffrey degrees of freedom
ci              # Bell-McCaffrey confidence interval
##### Save point estimates, SE, and t-stat
estimates.helpful <- cbind(point.estimate,bm$se,point.estimate/bm$se)[2:4,]
estimates.helpful # These results are reported in Appendix B in Table 4

### Create Figure 5
par(mfrow=c(1,1))
var.names <- c("Shaming", "Norms", "Norms and Shaming")
par(
  oma = c(0,0,0,0), # Since it is a single plot, I set the outer margins to zero.
  mar = c(5,4,5,4) # margins adjusted to reflect changing the locations of the labels
)
# create an empty plot for total customization
plot(NULL, # create empty plot
     ylim = c(-.3, .3), # set xlim by guessing
     xlim = c(1, .4+length(var.names)), # set ylim by the number of variables
     axes = F, xlab = NA, ylab = NA)
abline(h=0, lty=3)
est <- as.numeric(estimates.helpful[,1]) # conveniently store the estimates (minus the constant)
se <- as.numeric(estimates.helpful[,2]) # conveniently store the std. errors (minus the constant)
for (i in 1:length(est)) { # loop over a counter the length of the estimate vector
  lines(c(.2+i, .2+i), c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), lwd = 2, col="#B2CAF6")
  lines(c(.2+i, .2+i), c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), lwd = 4, col="#B2CAF6")
  points(.2+i, est[i],  pch = 19, cex = 2, col="white")
  points(.2+i, est[i], pch = 19, cex = 1, col="black")
  text(.2+i, est[i], var.names[i], xpd = T, cex = .8, pos = 3, offset=5) # add variable labels above the points
}
axis(side = 2) # add bottom axis
mtext(side = 2, "Estimated Treatment Effects", line = 2.5)

#########################
### Robustness Checks ###
#########################

### This section contains auxillary results described in the manuscript.

### Estimate logit and probit models mentioned in Footnote 9
#### Logit models
reply.logit <- glm(dat.sub$reply_prisoner ~ as.factor(dat.sub$sanctioning_no_norms) +
                     as.factor(dat.sub$norms_no_sanctioning) + as.factor(dat.sub$norms_sanctioning) +
                     as.factor(dat.sub$sender_two) + as.factor(dat.sub$wave) +
                     dat.sub$type + as.factor(dat.sub$ccode_host) + as.factor(dat.sub$separator) +
                     as.factor(dat.sub$separator_two) + as.factor(dat.sub$valediction), family=binomial(link='logit'))
summary(reply.logit)

helpful.logit <- glm(dat.sub$helpful_high ~ as.factor(dat.sub$sanctioning_no_norms) +
                     as.factor(dat.sub$norms_no_sanctioning) + as.factor(dat.sub$norms_sanctioning) +
                     as.factor(dat.sub$sender_two) + as.factor(dat.sub$wave) +
                     dat.sub$type + as.factor(dat.sub$ccode_host) + as.factor(dat.sub$separator) +
                     as.factor(dat.sub$separator_two) + as.factor(dat.sub$valediction), family=binomial(link='logit'))
summary(helpful.logit)

#### Probit models
reply.probit <- glm(dat.sub$reply_prisoner ~ as.factor(dat.sub$sanctioning_no_norms) +
                     as.factor(dat.sub$norms_no_sanctioning) + as.factor(dat.sub$norms_sanctioning) +
                     as.factor(dat.sub$sender_two) + as.factor(dat.sub$wave) +
                     dat.sub$type + as.factor(dat.sub$ccode_host) + as.factor(dat.sub$separator) +
                     as.factor(dat.sub$separator_two) + as.factor(dat.sub$valediction), family=binomial(link='probit'))
summary(reply.probit)

helpful.probit <- glm(dat.sub$helpful_high ~ as.factor(dat.sub$sanctioning_no_norms) +
                       as.factor(dat.sub$norms_no_sanctioning) + as.factor(dat.sub$norms_sanctioning) +
                       as.factor(dat.sub$sender_two) + as.factor(dat.sub$wave) +
                       dat.sub$type + as.factor(dat.sub$ccode_host) + as.factor(dat.sub$separator) +
                       as.factor(dat.sub$separator_two) + as.factor(dat.sub$valediction), family=binomial(link='probit'))
summary(helpful.probit)

#### Conduct permutation test mentioned on page 17; continue to find that the Norms cue is statistically significant
reply.fit.ri <- aovp(reply_prisoner ~ as.factor(sanctioning_no_norms) +
            as.factor(norms_no_sanctioning) + as.factor(norms_sanctioning), data=dat.sub)
summary(reply.fit.ri)

#### Conduct permutation test mentioned on page 19; continue to find that the Norms cue is statistically significant
helpful.fit.ri <- aovp(helpful_high ~ as.factor(sanctioning_no_norms) +
                      as.factor(norms_no_sanctioning) + as.factor(norms_sanctioning), data=dat.sub)
summary(helpful.fit.ri)

### Examine treatment effect heterogeneity
#### Check if treatment effects vary as a function of human rights scores
ols.fit.hrs <-  lm(dat.sub$reply_prisoner ~ as.factor(dat.sub$sanctioning_no_norms)*dat.sub$latentmean
               + as.factor(dat.sub$norms_no_sanctioning)*dat.sub$latentmean
               + as.factor(dat.sub$norms_sanctioning)*dat.sub$latentmean
               + as.factor(dat.sub$sender_two)
               + as.factor(dat.sub$wave)
               + dat.sub$type
               + as.factor(dat.sub$ccode_host)
               + as.factor(dat.sub$separator)
               + as.factor(dat.sub$separator_two)
               + as.factor(dat.sub$valediction), singular.ok = TRUE)
summary(ols.fit.hrs)

#### Check if treatment effects vary as a function of regime type
ols.fit.dem <-  lm(dat.sub$reply_prisoner ~ as.factor(dat.sub$sanctioning_no_norms)*dat.sub$democracy
               + as.factor(dat.sub$norms_no_sanctioning)*dat.sub$democracy
               + as.factor(dat.sub$norms_sanctioning)*dat.sub$democracy
               + as.factor(dat.sub$sender_two)
               + as.factor(dat.sub$wave)
               + dat.sub$type
               + as.factor(dat.sub$ccode_host)
               + as.factor(dat.sub$separator)
               + as.factor(dat.sub$separator_two)
               + as.factor(dat.sub$valediction), singular.ok = TRUE)
summary(ols.fit.dem)

#### Check if treatment effects vary as a function of UN region
ols.fit.reg <-  lm(dat.sub$reply_prisoner ~ as.factor(dat.sub$sanctioning_no_norms)*as.factor(dat.sub$un_region)
                   + as.factor(dat.sub$norms_no_sanctioning)*as.factor(dat.sub$un_region)
                   + as.factor(dat.sub$norms_sanctioning)*as.factor(dat.sub$un_region)
                   + as.factor(dat.sub$sender_two)
                   + as.factor(dat.sub$wave)
                   + dat.sub$type
                   + as.factor(dat.sub$ccode_host)
                   + as.factor(dat.sub$separator)
                   + as.factor(dat.sub$separator_two)
                   + as.factor(dat.sub$valediction), singular.ok = TRUE)
summary(ols.fit.reg)
